Vorgehensmodell
Business Understanding
Data Understandning
Modeling - Forecast der Ankuftszeiten
- Feature Engineering und Vorverarbeitung
  - One-Hot Encoding
  - Train-Test Split
- Modellauswahl
Evaluation & Deployment

Vorgehensmodell

Business Understanding

Für das Supply Chain Management der Daimler AG existiert ein Track und Trace System, das bestimmte Daten aufzeichnet. Ziel ist es aus den Daten nützliches Wissen für die Fachbereiche zu generieren, sodass die Lieferungen besser gesteuert werden können. Im ersten Schritt werden die Daten einer explorativen Datenanalyse unterzogen, um erste Erkenntnisse über die Lieferungen zu bekommen. Im Anschluss werden auf Top-Level Ebene einige Modelle entwickelt, welche die Ankuftszeiten vorhersagen. Dabei handelt es sich um erste Modelle, um die Einsatzmöglichkeiten zu evaluieren. Langfristig werden weiterhin diese Ziele verfolgt: -Forecast der Transportzeiten zwischen den einzelnen Stationen -Demand in den Produktionsstätten (entweder forecasten oder falls vorhanden mit aufnehmen) und anschließend mit aktuellen Lieferungen abgleichen, ob diese rechtzeitig da sind. -Klassifikation von Lieferungen in Risikogruppen

Data Understandning

Packages laden

library(tidyverse)
library(ggmap)
library(plotly)
library(dplyr)
library(lubridate)
library(ggridges)
library(ggplot2)
library(RColorBrewer)
library(bupaR)
library(scales)
library(mltools)
library(caret)
library(data.table)
library(mlr)
library(xgboost)
library(Metrics)
library(MLmetrics)
library(keras)
library(randomForest)
library(e1071)

Daten einlesen

raw_data <- read.csv("~/R 3.5 Files/Fallstudie Logistik/data_sample_case_study.csv", sep = ";")
raw_event_location <- read.csv("~/R 3.5 Files/Fallstudie Logistik/event_locations_case_study.csv", sep = ";")

Variablen untersuchen

#Variablen
colnames(raw_data)

##  [1] "ContainerNumber"  "ShipmentNumber"   "PoCreationDate"  
##  [4] "Material"         "ShipmentQuantity" "Controller"      
##  [7] "EventName"        "EventTime"        "EventMessageTime"
## [10] "EventLocation"    "VesselName"       "PortofDischarge" 
## [13] "Time2Arrival"

#Überblick der Daten
head(raw_data)

##   ContainerNumber ShipmentNumber   PoCreationDate    Material
## 1     HLXU8170197        4868820 19.02.2018 00:00 A2760101414
## 2     HLXU8170197        4868820 19.02.2018 00:00 A2760101414
## 3     HLXU8170197        4868820 19.02.2018 00:00 A2760101414
## 4     HLXU8170197        4868820 19.02.2018 00:00 A2760101414
## 5     HLXU8170197        4868820 19.02.2018 00:00 A2760101414
## 6     HLXU8170197        4868820 19.02.2018 00:00 A2760101414
##   ShipmentQuantity Controller       EventName        EventTime
## 1               72        T15    Loading Ship 07.04.2018 13:58
## 2               72        T15  Departure Ship 07.04.2018 19:05
## 3               72        T15  Unloading Ship 21.04.2018 06:19
## 4               72        T15 Departure Truck 21.04.2018 11:42
## 5               72        T15   Arrival Truck 08.05.2018 18:12
## 6               72        T15   Arrival Truck 08.05.2018 18:12
##   EventMessageTime            EventLocation  VesselName PortofDischarge
## 1 08.04.2018 03:19 urn:jaif:id:loc:26LNLRTM NYK ROMULUS                
## 2 09.04.2018 16:29 urn:jaif:id:loc:26LNLRTM NYK ROMULUS                
## 3 21.04.2018 11:59 urn:jaif:id:loc:26LUSCHS NYK ROMULUS                
## 4 21.04.2018 13:29 urn:jaif:id:loc:26LUSCHS NYK ROMULUS                
## 5 23.05.2018 22:04 urn:jaif:id:loc:26LUSCHS NYK ROMULUS                
## 6 23.05.2018 17:32 urn:jaif:id:loc:26LUSCHS NYK ROMULUS                
##   Time2Arrival
## 1          125
## 2          125
## 3          111
## 4          111
## 5           94
## 6           94

#Summary Shipment Quantity
summary(raw_data$ShipmentQuantity)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00   24.00   60.00   60.62   84.00  162.00

#Shipment Quantity Boxplot
plot_ly(y = raw_data$ShipmentQuantity, type = "box") %>%
  layout(title = "Boxplot Shipment Quantity",
         xaxis = list(title = ""),
         yaxis = list(title = ""))

Geocoding mit der Google API

Per Geocoding können die Koordinaten der Event Locations gezogen werden. Diese werden für eine spätere Visualisierung notwendig

register_google(key = "AIzaSyCwltxXK3VXCx0SDJxBPdfLMQWcEUOEvGk")
location <- as.character(raw_event_location$Bezeichnung)
cords <- geocode(location)
raw_event_location$lat <- cords$lat
raw_event_location$long <- cords$lon

Datenvorverarbeitung

#Left join, um die Bezeichnung sowie die Koordinaten in den gesamten Dataframe zu bekommen
data <- left_join(raw_data, raw_event_location, "EventLocation")

#Neue Spalte, die anzeigt wenn ein Teil angekommen ist
data$ShipmentArrived <- if_else(data$EventName == "Goods Receipt Dock" & data$Bezeichnung == "Werk Tuscaloosa", TRUE, FALSE)

#Umwandeln der Spaltentypen in Standardformat Zeit
data$EventTime <- as.POSIXct(as.character(data$EventTime),format = "%d.%m.%Y  %H:%M")
data$EventMessageTime <- as.POSIXct(as.character(data$EventMessageTime),format = "%d.%m.%Y  %H:%M")
data$PoCreationDate <- as.POSIXct(as.character(data$PoCreationDate),format = "%d.%m.%Y  %H:%M")
data$PoCreationDate <- as.Date(data$PoCreationDate)

#Als Numeric speichern
data$ShipmentNumber <- as.numeric(data$ShipmentNumber)

#Monat, Jahr, Wochentag und Quartal hinzufügen
data$monat <- month(data$PoCreationDate)
data$jahr <- year(data$PoCreationDate)
data$wochentag <- weekdays(data$PoCreationDate)
data$quater <- lubridate::quarter(data$PoCreationDate, with_year = TRUE, fiscal_start = 1)

#Daten 18.02.2018 bis 23.04.2019
summary(data$PoCreationDate)

##         Min.      1st Qu.       Median         Mean      3rd Qu. 
## "2018-02-18" "2018-06-25" "2018-09-25" "2018-09-24" "2019-01-08" 
##         Max. 
## "2019-04-23"

#Daten 05.04.2018 bis 24.05.2019
summary(data$EventTime)

##                  Min.               1st Qu.                Median 
## "2018-04-05 14:01:00" "2018-08-19 11:51:45" "2018-11-20 13:16:00" 
##                  Mean               3rd Qu.                  Max. 
## "2018-11-17 20:47:29" "2019-03-03 16:09:15" "2019-05-24 22:10:00" 
##                  NA's 
##                   "4"

#4 Einträge von EventTime haben ein NA
#Vorgehen: Da EventTime und EventTimeSend meistens ziemlich gleich sind werden die NA's damit ersetzt
data$EventTime <- if_else(is.na(data$EventTime) == TRUE, data$EventMessageTime, data$EventTime)

#Testen, welche Spalten weitere NA's haben
colnames(data)[colSums(is.na(data)) > 0]

## [1] "Bezeichnung" "lat"         "long"

#Filtern auf NA's
bez <- data %>%
  filter(is.na(Bezeichnung) == TRUE)

#Es ist keine Event Location vorhanden -> deswegen gibt es kein Mapping

#Für die genauere Analyse wird eine Sendungsnummer der betroffenen Events exploriert
bsp <- data %>%
  filter(ShipmentNumber == 5120761)

#Es sieht so aus als wäre diese Zeile ein Duplikat in dem die Event Location nicht richtig aufgezeichnet wurde
#Aus diesem Grund werden nun alle Duplikate entfernt, basierend auf der ContainerNummer, der ShipmentNumber, des ProCreationDate, des EventName, des Material, des Controller sowie der EventTime
data <- data[!duplicated(data[1:8]),]

#Testen, ob noch NA drin sind
colnames(data)[colSums(is.na(data)) > 0]

## [1] "Bezeichnung" "lat"         "long"

#Es sind noch 53 NA's enthalten -> wieder überprüfen
#Filtern auf NA's
bez <- data %>%
  filter(is.na(Bezeichnung) == TRUE)
bsp <- data %>%
  filter(ShipmentNumber == 5281788)

#Das sind alles Event mit dem Name ContainerClosed -> Das sind alle, die auch ein NA haben. Dementsprechend kann keine eindeutige Zuordnung des Events stattfinden.
#Dies sollte mit dem zuständigen Fachbereich geklärt werden!
#Für jetzt werden diese Events als eigenstädnige Event Location behandelt. Es wird sicher später zeigen, wo die Warenströme nach diesem Event hingehen
data$Bezeichnung <- if_else(data$EventName == "Container Closed", "Container Dummy", as.character(data$Bezeichnung))

#Die restlichen NA's in Latitude und Longitude können erstmal so gelassen werden, da keine Adresse vorhanden ist

#Diese dataframes wieder löschen, um RAM zu sparen
remove(bez)
remove(bsp)

Plot: Monatliche Quantity

Dieser Plot zeigt die monatlichen Mengen an, die transportiert werden getrennt nach Jahren. Es ist zu erkennen, dass 2019 bisher weniger transportiert wurde als im Vorjahr. Das Maximum lag im Jahr 2018 im April.

#Die Daten auf Basis der Sendungen
df_sendungen <- data[!duplicated(data[2]),]

#Nur auf monatlicher Ebene
df_sendungen <- df_sendungen %>%
  mutate(Date = format(as.Date(df_sendungen$PoCreationDate), "%Y-%m"))

#Anzahl an Jahre -> durch dieses Vorgehen wird sichergestellt, dass dieser Plot für eine beliebige Anzahl an Jahren richtig angezeigt wird
Years = as.data.frame(unique(year(df_sendungen$PoCreationDate)))
colnames(Years) <- c("Jahre")

#Plot initialisieren
p <- plot_ly(type = 'scatter', mode = 'lines', line = list(shape = "spline")) %>%
  layout(title = "Monatliche Quantity in Einheiten",
         xaxis = list(title = ""),
         yaxis = list(title = ""))

#Eigenen Farbvektor erstellen
color_vec <- c('rgb(8,48,107)','rgb(205, 12, 24)','rgb(0, 0, 0)','rgb(128, 128, 128)','rgb(0, 255, 255)', 'rgb(218, 165, 32)','rgb(34, 139, 34)')

#Für jedes Jahr wird ein Trace erstellt
i = 1
while (i <= nrow(Years)) {
  
  #Auf Jahr i filtern
  df_neu <- df_sendungen %>%
    filter(jahr == Years$Jahr[i])
  
  #Quantity auf monatlicher Ebene aggregieren
  eda <- aggregate(df_neu$ShipmentQuantity, by=list(ds=df_neu$Date), FUN=sum)
  eda$ds <- paste(eda$ds, "01", sep = "-")
  
  eda$x <- as.numeric(eda$x)
  eda$ds <- as.Date(eda$ds)
  eda$month <- month(eda$ds)
  name <- unique(year(eda$ds))
  
  #zum Plot hinzufügen
  p <- p %>%
    add_trace(data = eda, x = ~month, y = ~x , name = paste("Quantity ",name), line = list(color = color_vec[i]))  
  
  i = i+1
}

p

Plot: Dauer von jedem Hafen bis nach Tuscaloosa

Dieser Plot zeigt die Dauer von jedem Hafen bis nach Tuscaloosa. Hier ist wohl angemerkt, dass die Daten nicht auf der Eventebene betrachtet werden. Das bedeutet, dass unterschiedliche Events die Daten verzerren

#Die gleichen Werte bei unterschiedlichen Bestellungen entfernt
df_p1 <- data[!duplicated(data[4:11]),]
  
df_p1 <- df_p1 %>%
  select(ShipmentNumber, Bezeichnung, Time2Arrival)

#Hier wieder die Duplikate entfernen -> Sonst verzerrt
df_p1 <- df_p1[!duplicated(df_p1[1:3]),]


plot_ly(df_p1, y = ~Time2Arrival, color = ~Bezeichnung, type = "box")

Plot: Dauer von jedem Hafen inklusive Events bis nach Tuscaloosa

Dieser Plot zeigt die Dauer von jedem Hafen bis nach Tuscaloosa. Hier ist die Detailebene genauer, da die Darstellung auf der Hafen und Eventebene liegt. Es ist deutlich zu erkennen, dass bestimmte Häfen sehr hohe Ausreißer haben (z.B. Bremerhafen oder Charleston)

#Die gleichen Werte bei unterschiedlichen Bestellungen entfernt
df_p1 <- data[!duplicated(data[4:11]),]

#Hafen und Events verbinden
df_p1$HafenEvent <- paste(df_p1$Bezeichnung, df_p1$EventName,  sep = "-")
  
df_p1 <- df_p1 %>%
  select(ShipmentNumber, HafenEvent, Time2Arrival)

#Hier wieder die Duplikate entfernen -> Sonst verzerrt
df_p1 <- df_p1[!duplicated(df_p1[1:3]),]

plot_ly(df_p1, y = ~Time2Arrival, color = ~HafenEvent, type = "box") %>%
  layout(showlegend = FALSE)

Plot: Häufigkeit verwendeter Schiffe

Dieser Plot zeigt die Häufigkeit der verwendeten Schiffe. Die ersten drei Schiffe werden ungefähr gleich häufig verwendet. Das spätere Modell wird zeigen, welchen Einfluss die Verwendung der Schiffe auf die Ankunftszeit hat.

#Häufigkeit der Schiffe
df_p2 <- count(df_sendungen, vars = VesselName)

#Absteigend ordnen und als Factor speichern
df_p2 <- df_p2[order(df_p2$n),] 
df_p2$vars <- factor(df_p2$vars, levels = c(as.character(df_p2$vars)))

#Plottet die Häufigkeit der Schiffe
plot_ly(df_p2, x = ~n, y = ~vars, type = 'bar', 
        marker = list(color = 'rgb(158,202,225)',
                      line = list(color = 'rgb(8,48,107)', width = 1.5))) %>%
  layout(title = "Häufigkeit der Schiffe",
         xaxis = list(title = ""),
         yaxis = list(title = ""))

Plot: Materialbedarf nach Quartalen

Dieser Plot zeigt den Materialbedarf in den einzelnen Quartalen an. Es ist zu erkennen, dass der Materialbedarf im 1. Quartal 2019 leicht über dem des Vorjahres liegt. Mit mehr Daten könnten hier Zusammenhänge besser interpretiert werden.

#Daten nach Quartal gruppieren und die Mengen summieren
df_p3 <- df_sendungen %>% 
  group_by(quater) %>% 
  summarise(Quantity = sum(ShipmentQuantity))

#Achsen richtig benennen und den Durchschnitt als Vergleichsgrundlage berechnen
df_p3$QuaterName <- paste("Quartal", df_p3$quater, sep = " ")
df_p3$avg <- mean(df_p3$Quantity)

p <- plot_ly(df_p3) %>%
  add_trace(x = ~QuaterName, y = ~Quantity, type = 'bar', name = 'Materialbedarf',
            marker = list(color = 'rgb(158,202,225)',
                          line = list(color = 'rgb(8,48,107)', width = 1.5))) %>%
  add_trace(x = ~QuaterName, y = ~avg, type = 'scatter', mode = 'lines', name = 'Durchschnitt',
            line = list(color = 'rgb(205, 12, 24)',width = 0.7),
            hoverinfo = "text") %>%
  layout(title = "Materialbedarf nach Quartalen",
         bargap = 0.4,
         xaxis = list(title = ""),
         yaxis = list(title = ""),
         legend = list(orientation = 'h'))

p

Plot: Density nach Materialart

Dieser Plot zeigt die Verteilung der monatlichen Quantity Mengen nach Materialart an. Insbesondere A1662708602 hat starke Ausreißer nach oben.

theme_set(theme_ridges())

#Nach Date und Materialart gruppieren und summieren
df_p4 <- df_sendungen %>% 
  group_by(Date, Material) %>% 
  summarise(Quantity = sum(ShipmentQuantity))

 
#Die schönen Paletten haben nur 9 Farben -> 1 muss künstlich hinzugefügt werden
mycolors <- brewer.pal(n = 9, name = "Blues") 
mycolors <- c(mycolors, "#001642")

ggplot(df_p4, aes(x = Quantity, y = Material)) +
  geom_density_ridges(aes(fill = Material)) +
  scale_fill_manual(values = mycolors)+
  labs(title = 'Verteilung der monatlichen Quantity Mengen nach Materialart')

Supply Chain auf der Hafenebene

Plot: Supply Chain Map - Häufigkeit

Dieser Plot zeigt die Struktur der zugrundeliegenden Supply Chain an. Dafür werden Techniken des Process Mining mit dem Bupar Package eingesetzt. Hierbei fungiert die ShipmentNumber als CaseID, der Hafen als Event und das Schiff als Resource. Hier ist zu erkennen wie die Vernetzung zwischen den Häfen ist.

df_pm <- data

#Die Aktivität wird definiert als Hafenbezeichnung
df_pm$Activity <- df_pm$Bezeichnung

#Nur diese Spalten selektieren
df_pm <- df_pm %>%
  select(ShipmentNumber, Activity, VesselName, EventTime)

#Da die Daten auf der granularen Ebene der Bestellungen ist, müssen die Duplikate entfernt werden -> wir wollen nur die Sendungen verfolgen
df_pm <- unique(df_pm)

#Künstliche Activity Instance hinzufügen
df_pm <- df_pm %>%
  mutate(activity_instance = 1:nrow(.))

#Künstlichen LifeCycle hinzufügen, da keine Informationen vorhanden sind für die Länge der Events
df_pm$LifeCycle <- "End"

#Duplikate entfernen
df_pm <- df_pm[!duplicated(df_pm[1:2]),]

#Event Log kreieren
event_logs <- df_pm %>% 
  eventlog(
    case_id = "ShipmentNumber",
    activity_id = "Activity",
    activity_instance_id = "activity_instance",
    lifecycle_id = "LifeCycle",
    timestamp = "EventTime",
    resource_id = "VesselName"
  )


#Prozess Map plotten
event_logs %>%
  filter_activity_frequency(percentage = .99) %>% 
  filter_trace_frequency(percentage = .99) %>%    
  process_map(type = frequency("absolute"))

Plot: Supply Chain Map - Durchschnittliche Dauer

event_logs %>%
      filter_activity_frequency(percentage = .99) %>% 
      filter_trace_frequency(percentage = .99) %>%    
      process_map(processmapR::performance(mean, "days"))

Plot: Vorgänger Nachgänger Beziehung der Häfen

event_logs %>%
  precedence_matrix(type = "relative_consequent") %>%
  plot

Plot: Durchschnittliche Zeit nach Schiffen

event_logs %>%
  idle_time("resource", units = "days") %>%
  plot()

Plot: Weltkarte der Transportwege

Diese Weltkarte zeigt die Transportwege der Supply Chain. Die Größte der Punkte steht für die ankommende Menge an diesem Ort.

#Die Vorgänger Nachgänger Beziehungen extrahieren, um die Transportwege zu visualisieren
pred_matrix <- event_logs %>%
  precedence_matrix(type = "relative_consequent")

#Diese Verbindungen haben keinen definierten Ort
df_map <- pred_matrix %>%
  filter(consequent != "End") %>%
  filter(antecedent != "Start")


#Daten joinen, um lat und long zu bekommen
df_map <- left_join(df_map, raw_event_location, by = c("antecedent" = "Bezeichnung"))
df_map <- df_map %>%
  select(antecedent, consequent, rel_consequent, n, lat, long)
colnames(df_map) <- c("antecedent", "consequent", "rel_consequent", "n", "lat_antecedent", "long_antecedent")

df_map <- left_join(df_map, raw_event_location, by = c("consequent" = "Bezeichnung"))
df_map <- df_map %>%
  select(antecedent, consequent, rel_consequent,n, lat_antecedent, long_antecedent, lat, long)
colnames(df_map) <- c("antecedent", "consequent", "rel_consequent","n", "lat_antecedent", "long_antecedent", "lat_consequent", "long_consequent")

#Zur Ansicht die Daten skalieren
df_map$nscale <- rescale(df_map$n, to=c(1,5))

#Size nach Anzahl der ankommenden Lieferungen und der ausgehenden Lieferungen -> Filtermöglichkeiten für später
ag <- aggregate(list(AnkommendeMenge=df_map$n), by=list(Bezeichnung=df_map$consequent), FUN=sum)
raw_event_location <- left_join(raw_event_location, ag, "Bezeichnung")

ag <- aggregate(list(AusgehendeMenge=df_map$n), by=list(Bezeichnung=df_map$antecedent), FUN=sum)
raw_event_location <- left_join(raw_event_location, ag, "Bezeichnung")

#NA mit 0 ersetzen
raw_event_location$AnkommendeMenge[is.na(raw_event_location$AnkommendeMenge)] <- 0

#Die Mengen skalieren, da sonst zu groß
raw_event_location$ScaledAnkommendeMenge <- rescale(raw_event_location$AnkommendeMenge, to=c(0.25,2))
raw_event_location$ScaledAusgehendeMenge <- rescale(raw_event_location$AusgehendeMenge, to=c(0.25,2))


#Fahrt benennen, um den hover Text anzupassen
df_map$FahrtName <- paste(df_map$antecedent,df_map$consequent, sep = " - ")

geo <- list(
  scope = 'world',
  projection = list(type = 'azimuthal equal area'),
  showland = TRUE,
  landcolor = toRGB("gray95"),
  countrycolor = toRGB("gray80")
)

p <- plot_geo(locationmode = 'USA-states')

for(i in 1:nrow(raw_event_location)){
  df_el <- as.data.frame(raw_event_location[i,])
  p <- p %>%
    add_markers(
    data = df_el, x = ~long, y = ~lat, text = ~Bezeichnung,
    size = ~ScaledAnkommendeMenge, hoverinfo = "text", alpha = 0.5
  )
}
p <- p %>%
  add_segments(
    data = df_map,
    x = ~long_antecedent, xend = ~long_consequent,
    y = ~lat_antecedent, yend = ~lat_consequent, text = ~FahrtName,
    alpha = 0.3, size = ~nscale, hoverinfo = "text"
  ) %>%
  layout(
    title = 'Supply Chain Routen - Mercedes Benz',
    geo = geo, showlegend = FALSE
  )

p

Modeling - Forecast der Ankuftszeiten

Feature Engineering und Vorverarbeitung

Nach reiflicher Überlegung wurden 7 Features ausgewählt, die einen potenziellen Einfluss auf die Ankuftszeit haben. Diese sind: -Material -EventName -EventLocation bzw. Bezeichnung -VesselName -Monat -Jahr -Wochentag

#Zuerst muss der Datensatz auf die ShipmentNumber reduziert werden -> Auf dieser Basis werden die Bestellungen rausgefiltert.
#Ansonsten hat man doppelte Daten drin
df <- data %>%
  select(ShipmentNumber, Material, EventName, Bezeichnung, VesselName, monat, jahr, wochentag, ShipmentArrived, Time2Arrival)
df <- df[!duplicated(df[1:9]),]

#Herausfiltern von ShipmentArrived = TRUE. Macht keinen Sinn das forezucasten
df <- df %>%
  filter(ShipmentArrived == FALSE)

#ShipmentArrived und ShipmentNumber droppen
df <- df %>%
  select(-ShipmentArrived, -ShipmentNumber)

One-Hot Encoding

Mit One-Hot Encoding werden alle Categorial Variablen umgewandelt in numerische. Damit können die Machine Learning Verfahren umgehen.

colnames(df)

## [1] "Material"     "EventName"    "Bezeichnung"  "VesselName"  
## [5] "monat"        "jahr"         "wochentag"    "Time2Arrival"

#Als Factor speichern
df$Material <- as.factor(df$Material)
df$EventName <- as.factor(df$EventName)
df$Bezeichnung <- as.factor(df$Bezeichnung)
df$VesselName <- as.factor(df$VesselName)
df$monat <- as.factor(df$monat)
df$jahr <- as.factor(df$jahr)
df$wochentag <- as.factor(df$wochentag)

#One-Hot Encoding für alle Spalten
df_tb <- setDT(df)
df_tb <- one_hot(df_tb, cols = "Material")
df_tb <- one_hot(df_tb, cols = "EventName")
df_tb <- one_hot(df_tb, cols = "Bezeichnung")
df_tb <- one_hot(df_tb, cols = "VesselName")
df_tb <- one_hot(df_tb, cols = "monat")
df_tb <- one_hot(df_tb, cols = "jahr")
df_tb <- one_hot(df_tb, cols = "wochentag")

#Wir haben jetzt 81 Features
df <- as.data.frame(df_tb)
colnames(df)

##  [1] "Material_A1662707202"                   
##  [2] "Material_A1662708602"                   
##  [3] "Material_A1662708702"                   
##  [4] "Material_A2132704801"                   
##  [5] "Material_A2560107300"                   
##  [6] "Material_A2572700000"                   
##  [7] "Material_A2760101414"                   
##  [8] "Material_A2760101614"                   
##  [9] "Material_A2760106714"                   
## [10] "Material_A6540107406"                   
## [11] "EventName_Arrival Ship"                 
## [12] "EventName_Arrival Truck"                
## [13] "EventName_Container Closed"             
## [14] "EventName_Departure Ship"               
## [15] "EventName_Departure Truck"              
## [16] "EventName_Goods Receipt Dock"           
## [17] "EventName_Loading Ship"                 
## [18] "EventName_Loading Truck"                
## [19] "EventName_New Scheduling"               
## [20] "EventName_Registration Yard"            
## [21] "EventName_Unloading Ship"               
## [22] "Bezeichnung_Bahnterminal Savannah"      
## [23] "Bezeichnung_Consolidation Center Speyer"
## [24] "Bezeichnung_Container Dummy"            
## [25] "Bezeichnung_Hafen Aguascalientes"       
## [26] "Bezeichnung_Hafen Altamira"             
## [27] "Bezeichnung_Hafen Antwerpen"            
## [28] "Bezeichnung_Hafen Bremerhaven"          
## [29] "Bezeichnung_Hafen Charleston"           
## [30] "Bezeichnung_Hafen Hamburg"              
## [31] "Bezeichnung_Hafen Rotterdam"            
## [32] "Bezeichnung_Werk Tuscaloosa"            
## [33] "VesselName_BREVIK BRIDGE"               
## [34] "VesselName_CHARLESTON EXPRESS"          
## [35] "VesselName_EASLINE DALIAN"              
## [36] "VesselName_MAINE TRADER"                
## [37] "VesselName_MOL EMISSARY"                
## [38] "VesselName_MOL EMPIRE"                  
## [39] "VesselName_MOL EXPERIENCE"              
## [40] "VesselName_MONACO MAERSK"               
## [41] "VesselName_MUNICH MAERSK"               
## [42] "VesselName_NYK DAEDALUS"                
## [43] "VesselName_NYK DENEB"                   
## [44] "VesselName_NYK METEOR"                  
## [45] "VesselName_NYK NEBULA"                  
## [46] "VesselName_NYK REMUS"                   
## [47] "VesselName_NYK RIGEL"                   
## [48] "VesselName_NYK ROMULUS"                 
## [49] "VesselName_NYK RUMINA"                  
## [50] "VesselName_PHILADELPHIA EXPRESS"        
## [51] "VesselName_ROTTERDAM EXPRESS"           
## [52] "VesselName_SHANGHAI TRADER"             
## [53] "VesselName_ST LOUIS EXPRESS"            
## [54] "VesselName_TRF PARTICI"                 
## [55] "VesselName_VECCHIO BRIDGE"              
## [56] "VesselName_VENICE BRIDGE"               
## [57] "VesselName_WASHINGTON EXPRESS"          
## [58] "VesselName_YM UNITY"                    
## [59] "VesselName_YORKTOWN EXPRESS"            
## [60] "monat_1"                                
## [61] "monat_2"                                
## [62] "monat_3"                                
## [63] "monat_4"                                
## [64] "monat_5"                                
## [65] "monat_6"                                
## [66] "monat_7"                                
## [67] "monat_8"                                
## [68] "monat_9"                                
## [69] "monat_10"                               
## [70] "monat_11"                               
## [71] "monat_12"                               
## [72] "jahr_2018"                              
## [73] "jahr_2019"                              
## [74] "wochentag_Dienstag"                     
## [75] "wochentag_Donnerstag"                   
## [76] "wochentag_Freitag"                      
## [77] "wochentag_Mittwoch"                     
## [78] "wochentag_Montag"                       
## [79] "wochentag_Samstag"                      
## [80] "wochentag_Sonntag"                      
## [81] "Time2Arrival"

remove(df_tb)

Train-Test Split

Um das Modell anschließend evaluieren zu können, wird der Datensatz unterteilt in einen Trainingsset und Testset (80/20)

#Definition der Sample size
smp_size <- floor(0.8 * nrow(df))

#Set Seed
set.seed(100)
train_ind <- sample(seq_len(nrow(df)), size = smp_size)

train <- df[train_ind, ]
test <- df[-train_ind, ]

Modellauswahl

Da relativ wenige Daten vorhanden sind, werden voraussichtlich neuronale Netze nicht die besten Ergebnisse liefern. Sie werden hier trotzdem mit aufgeführt, um für die Skalierbarkeit bereits das Modell zu haben.

XGBoost - Modell

#Als Datatable speichern
train_x <- setDT(train)
test_x <- setDT(test)

#Targetvariablen speichern
labels_train <- train_x$Time2Arrival
labels_test <- test_x$Time2Arrival

#Targetvariablen entfernen, sodass nur Features da sind
train_x <- subset(train_x, select = - Time2Arrival)
test_x <- subset(test_x, select = - Time2Arrival)

#Als Matrix speichern
train_x <- data.matrix(train_x, rownames.force = NA)
test_x <- data.matrix(test_x, rownames.force = NA)

#Als XGBoost MAtrix speichern
dtrain <- xgb.DMatrix(data = train_x,label = labels_train) 
dtest <- xgb.DMatrix(data = test_x,label=labels_test)

#Für die Evaluation während des Trainings
watchlist <- list(train=dtrain, test=dtest)
set.seed(101)

st <- xgb.train(data=dtrain, booster = "gbtree", max.depth=7, lambda = 2, alpha = 2, colsample_bytree = 0.5, min_child_weight = 0.2, early_stopping_rounds= 200, nthread = 2, gamma = 0, eta= 0.01, nrounds=10000, watchlist=watchlist, eval.metric = "rmse", objective = "reg:linear")

#Auf die Testdaten predicten
xgbpred <- predict (st,dtest)

#Für späteres Stacking
xgbpred_train <- predict(st, dtrain)

#xgb.save(st, "xgb.model")

XGBoost - Validation

Das vorher erstellte Modell wird nun anhand der Testdaten validiert.

#Validation Kennzahlen
rmse(labels_test, round(xgbpred,0))

## [1] 17.62393

MAE(labels_test, round(xgbpred,0))

## [1] 9.39945

#Feautre Importance
mat <- xgb.importance(feature_names = colnames(train_x),model = st)
xgb.plot.importance (importance_matrix = mat[1:30])

#Ausschnitt des Tree's
model <- xgb.dump(st, with.stats = T)
model[1:30]

##  [1] "booster[0]"                                                    
##  [2] "0:[f19<0.5] yes=1,no=2,missing=1,gain=1348963,cover=20366"     
##  [3] "1:[f28<0.5] yes=3,no=4,missing=3,gain=848633,cover=17549"      
##  [4] "3:[f71<0.5] yes=7,no=8,missing=7,gain=290917,cover=8785"       
##  [5] "7:[f21<0.5] yes=15,no=16,missing=15,gain=33451.75,cover=1980"  
##  [6] "15:[f5<0.5] yes=31,no=32,missing=31,gain=6525.875,cover=1784"  
##  [7] "31:[f3<0.5] yes=57,no=58,missing=57,gain=2595.375,cover=1537"  
##  [8] "57:[f33<0.5] yes=93,no=94,missing=93,gain=1005.375,cover=1321" 
##  [9] "93:leaf=0.280295104,cover=1184"                                
## [10] "94:leaf=0.327302158,cover=137"                                 
## [11] "58:leaf=0.237935767,cover=216"                                 
## [12] "32:leaf=0.217811242,cover=247"                                 
## [13] "16:[f45<0.5] yes=33,no=34,missing=33,gain=1063.27344,cover=196"
## [14] "33:[f5<0.5] yes=59,no=60,missing=59,gain=2070.19922,cover=183" 
## [15] "59:[f46<0.5] yes=95,no=96,missing=95,gain=173.857422,cover=177"
## [16] "95:leaf=0.133896098,cover=152"                                 
## [17] "96:leaf=0.0850000009,cover=25"                                 
## [18] "60:leaf=0.308750004,cover=6"                                   
## [19] "34:leaf=0.0350000001,cover=13"                                 
## [20] "8:[f21<0.5] yes=17,no=18,missing=17,gain=158137,cover=6805"    
## [21] "17:[f30<0.5] yes=35,no=36,missing=35,gain=124457,cover=6551"   
## [22] "35:[f5<0.5] yes=61,no=62,missing=61,gain=31794,cover=6422"     
## [23] "61:[f3<0.5] yes=97,no=98,missing=97,gain=31640,cover=5675"     
## [24] "97:leaf=0.41574043,cover=5049"                                 
## [25] "98:leaf=0.336449057,cover=626"                                 
## [26] "62:[f29<0.5] yes=99,no=100,missing=99,gain=5413.625,cover=747" 
## [27] "99:leaf=0.344053268,cover=674"                                 
## [28] "100:leaf=0.236733332,cover=73"                                 
## [29] "36:[f60<0.5] yes=63,no=64,missing=63,gain=36895.6875,cover=129"
## [30] "63:leaf=0.599607825,cover=100"

Neuronales Netz - Modell

train_labels <- labels_train
train_data <- as.data.frame(train_x)

test_labels <- labels_test
test_data <- as.data.frame(test_x)

#Dataframe als numeric speichern
train_data  <- as.data.frame(lapply(train_data, as.numeric))
test_data  <- as.data.frame(lapply(test_data, as.numeric))

#Check ob wirklich alle NA's raus sind -> Neuronale Netze können damit nicht umgehen
colnames(train_data)[colSums(is.na(train_data)) > 0]

## character(0)

colnames(test_data)[colSums(is.na(test_data)) > 0]

## character(0)

#Als Matrix speichern
train_data <- as.matrix(train_data)
test_data <- as.matrix(test_data)

#Neuronales Netz mit einem Input Layer, ein Hidden Layer und ein Output Layer
build_model <- function() {
  
  model <- keras_model_sequential() %>%
    layer_dense(units = 64, activation = "relu", input_shape = dim(train_data)[2]) %>%
    layer_dense(units = 32, activation = "relu") %>%
    layer_dense(units = 1)
  
  model %>% compile(
    loss = 'mean_squared_error',
    optimizer = optimizer_rmsprop(),
    metrics = list("mean_absolute_error")
  )
  
  model
}

model <- build_model()
model %>% summary()

## ___________________________________________________________________________
## Layer (type)                     Output Shape                  Param #     
## ===========================================================================
## dense (Dense)                    (None, 64)                    5184        
## ___________________________________________________________________________
## dense_1 (Dense)                  (None, 32)                    2080        
## ___________________________________________________________________________
## dense_2 (Dense)                  (None, 1)                     33          
## ===========================================================================
## Total params: 7,297
## Trainable params: 7,297
## Non-trainable params: 0
## ___________________________________________________________________________

#Funktion damit das Modell die Scores während des Trainings ausgibt
print_dot_callback <- callback_lambda(
  on_epoch_end = function(epoch, logs) {
    if (epoch %% 80 == 0) cat("\n")
    cat(".")
  }
)    

#Rundenanzahl
epochs <- 500

#Early Stopping wenn sich der Score nicht verbessert
early_stop <- callback_early_stopping(monitor = "val_loss", patience = 20)

model <- build_model()
history <- model %>% fit(
  train_data,
  train_labels,
  epochs = epochs,
  validation_split = 0.2,
  verbose = 0,
  batch_size = 32,
  callbacks = list(early_stop, print_dot_callback)
)

## 
## ................................................................................
## ........................

test_predictions <- model %>% predict(test_data)
train_predictions_nn <- model %>% predict(train_data)

Neuronales Netz - Validation